Лекция 2. Решения уравнения фильтрации

версия 0.2 от 17.09.2025

Хабибуллин Ринат

# импортируем библиотеки, которые могут пригодиться для проведения расчетов
import numpy as np
import matplotlib.pyplot as plt
import scipy.special as spec
import anaflow
import sympy as sp
import mpmath as mp

1 Уравнение фильтрации и его решения в пространстве Лапласа

Рассматривается уравнение фильтрации в безразмерных переменных

\[ \frac{\partial p_D}{ \partial t_D} = \frac{1}{r_D}\left[ \frac{ \partial{}}{ \partial{r_D} }\left( r_D \dfrac{\partial p_D}{ \partial r_D} \right) \right] \tag{1.1}\]

где введены следующие безразмерные переменные

• \(r_D\) - безразмерное расстояние от центра скважины
• \(t_D\) - безразмерное время
• \(p_D\) - безразмерное давление
• \(q_D\) - безразмерный дебит

Соответствующие определения безразмерных переменных

\[ r_D = \frac{r}{r_w} \tag{1.2}\]

\[ t_D = \frac{0.00036 kt}{\phi \mu c_t r_w^2} \tag{1.3}\] \[ p_D = \frac{kh}{ 18.42 q_s B \mu} \left( p_i - p \right) \tag{1.4}\]

где в свою очередь

• \(q_s\) - дебит скважины на поверхности, приведенный к нормальным условиям, м3/сут
• \(\phi\) - пористость, доли единиц
• \(\mu\) - вязкость нефти в пласте, сП
• \(B\) - объемный коэффициент нефти, м3/м3
• \(p_i\) - начальное давление в пласте, атм
• \(p\) - давление забойное, атм
• \(c_t\) - общая сжимаемость системы в пласте, 1/атм
• \(k\) - проницаемость, мД
• \(t\) - время, час
• \(r\) - расстояние от центра скважины, м
• \(r_w\) - радиус скважины, м

1.1 Общее решение уравнения фильтрации в пространстве Лапласа

Решение такого уравнение может быть получено с использованием преобразования Лапласа.

\[ L \left [ f(t) \right] = \tilde{f}(u) = \int_{0}^{\infty}f(t)e^{-ut}dt \]

где \(u\) параметр пространства Лапласа соответствующий времени.

Тогда уравнение в пространстве Лапласа преобразуется к виду:

\[ u \tilde{p}_D = \dfrac{1}{r_D} \left[\dfrac{d}{d r_D} \left(r_D \dfrac{d{\tilde{p}_D}}{d r_D} \right) \right] \tag{1.5}\]

Уравнение (1.5) известно как модифицированное уравнение Бесселя. Общее решение этого уравнения можно записать в виде

\[ \tilde{p}_D(u, r_D) = A(u) K_0(r_D \sqrt u) + B(u) I_0(r_D \sqrt u) \tag{1.6}\]

где

• \(u\) - переменная пространства Лапласа, соответствующая времени
• \(\tilde{p}_D(u, r_D)\) - изображение давления в пространстве Лапласа
• \(K_0, I_0\) - модифицированные функции Бесселя нулевого порядка (могут быть вычислены, например, с использованием реализации в библиотке scipy.special)
• \(A(u), B(u)\) - произвольные функции, которые могут быть определены при задании начальных и граничных условий

Для модифицированных функций Бесселя нулевого и первого порядка можно записать соотношения

\[ \dfrac{dI_0(u)}{du} = I_1(u) \tag{1.7}\]

\[ \dfrac{dK_0(u)}{du} =-K_1(u) \tag{1.8}\]

Для построения простых решений далее пригодятся некоторые свойства преобразования Лапласа

\[ \mathcal{L} \left [ a \right] = \dfrac{a}{u} \tag{1.9}\]

\[ \mathcal{L} \left [ \dfrac{df}{dt} \right] = u \tilde{f}(s) - f(t=0) \tag{1.10}\]

x = np.arange(1,5,0.1)

plt.rcParams["figure.figsize"] = (8,4)
fig, (ax1, ax2) = plt.subplots(1,2)
ax1.plot(x, spec.kn(0,x), label='$K_0$')
ax1.plot(x, spec.kn(1,x), label='$K_1$')
ax1.set_title("Bessel functions $K_0$, $K_1$")
ax1.legend()

ax2.plot(x, spec.iv(0,x), label='$I_0$')
ax2.plot(x, spec.iv(1,x), label='$I_1$')
ax2.set_title("Bessel functions $I_0$, $I_1$")
ax2.legend()

plt.tight_layout()
plt.show()

Рисунок 1.1: графики модифицированных функций Бесселя

1.2 Решение линейного стока - простейшее решение

Для построения частного решения необходимо задать начальные и граничные условия. Простейшее решение можно получить задав следующие начальные и граничные условия:

1. Однородное начальное давление \[ p_D(t_D=0, r_D) = 0 \]

2. Граничное условия на бесконечности

\[ \lim_{r_D \to \infty} p_D(r_D, t_D) = 0 \tag{1.11}\]

в пространстве Лапласа (1.11) преобразуется в следующее

\[ \lim_{r_D \to \infty} \tilde{p}_D = 0 \tag{1.12}\]

3. Граничное условие на скважине

\[ \lim_{r_D \to r_{wD}} \left[ r_D \dfrac{ \partial p_D(r_D, t_D)}{\partial r_D} \right] = -1 \tag{1.13}\]

в пространстве Лапласа (1.13) с учетом выражения (1.9) преобразуется к

\[ \lim_{r_D \to r_{wD}} \left[ r_D \dfrac{ d \tilde{p}_D}{d r_D} \right] = -\dfrac{1}{u} \tag{1.14}\]

Граничное условие на скважине (1.13) можно записать в размерном виде используя определение (1.4) и сравнить его с законом Дарси в радиальной форме \(q_{res} = \dfrac{kh}{18.42 \mu} r \dfrac{dP}{dr}\).

\[ \lim_{r \to r_{w}} \left[ r \dfrac{ \partial p(r, t)}{\partial r} \right] = 18.42 \dfrac{qB\mu}{kh} \]

Здесь надо учесть, что \(q_{res}\) это дебит из пласта при соответствующих далении и температуре, а \(q\) это дебит на поверхности приведенный к стандартным условиям. Они отличаются на объемный коэффициент \(B\).

Для построения частного решения необходимо исходя из приведенных условий подобрать значения функций \(A(u)\) и \(B(u)\) для общего решения (1.6).

\[ \tilde{p}_D(u, r_D) = A(u) K_0(r_D \sqrt u) + B(u) I_0(r_D \sqrt u) \]

Глядя на поведение функции \(I_0\) на бесконечности, понятно что для нашего случая надо положить \(B(u) = 0\).

Тогда решение имеет вид

\[ \tilde{p}_D = A(u) K_0(r_D \sqrt u) \tag{1.15}\]

Подставим выражение (1.15) в граничное условие (1.14)

\[ \lim_{r_D \to r_{wD}} \left[ r_D \dfrac{ d \left( A(u) K_0(r_D \sqrt u)\right)}{d r_D} \right] = -\dfrac{1}{u} \tag{1.16}\]

Преобразуем (1.16) учитывая (1.8)

\[ \lim_{r_D \to r_{wD}} \left[A(u) r_D\sqrt u K_1(r_D \sqrt u) \right] = \dfrac{1}{u} \]

Подставим вместо расстояния - предельное значение

\[ A(u) r_{wD} \sqrt u K_1(r_{wD} \sqrt u) = \dfrac{1}{u} \]

Наконец получим выражение для \(A(u)\)

\[ A(u) = \dfrac {1}{u r_{wD} \sqrt u K_1(r_{wD} \sqrt u)} \]

Тогда решение для произвольного радиуса будет иметь вид

\[ \tilde{p}_D(u) = \frac{K_0 \left( r_D \sqrt u \right)}{u r_{wD} \sqrt u K_1(r_{wD} \sqrt u)} \]

при \(r_{wD} = 1\) получим

\[ \tilde{p}_D(u) = \frac{K_0 \left( r_D \sqrt u \right)}{u \sqrt u K_1(\sqrt u)} \tag{1.17}\]

при \(r_{wD} = 0\), используя свойства модифицированной функции Бесселя \(\lim_{s \to 0} sK_1(s) = 1\) (смотри рисунок 1.1) получим \(A(u) = \dfrac{1}{u}\)

Решение для бесконечно малого радиуса скважины в пространстве Лапласа будет иметь вид

\[ \tilde{p}_D(u) = \frac{1}{u} K_0 \left( r_D \sqrt u \right) \tag{1.18}\]

где

\(K_0\), \(K_1\) - модифицированные функции Бесселя

Хотя решения в пространстве Лапласа относительно легко получить – обратная процедура получения решения в исходных координатах на основе решения в пространстве Лапласа оказывается сложнее. Аналитически это не всегда удается сделать, чаще эту процедуру проводят численно.

Для численного обратного преобразования Лапласа можно, например, использовать библиотеку mpmath Там же можно найти численную реализацию функций Бесселя. Но такой вариант расчета относительно медленный (смотрите пример 00_некоторые_технические_подробности.ipynb для подробностей). Дальнейшие расчеты будут основаны на scipy.special и anaflow.

1.3 Построение и исследование решения линейного стока с использованием обратного преобразования Лапласа

Численного обратное преобразование Лапласа удобный инструмент для построение решений уравнения фильтрации. Последующие примеры показывают как можно его использовать с python.

# пример функции реализующий расчет решения в пространстве Лапласа
def pd_lapl_line_source(u, rd=1.):
    """
    расчет решения линейного стока для безразмерного давления в пространстве Лапласа
    u - переменная пространства Лапласа
    rd- безразмерное расстояние от центра скважины
    """
    return np.divide(spec.kn(0, rd * u**0.5) , u)

# функция расчета безразмерного давления с использованием алгоритма Стефеста
# для численного обратного преобразования Лапласа
def pd_line_source_inv(td, rd=1.):
    """
    расчет решения линейного стока для безразмерного давления
    на основе численного обратного преобразования Лапласа (алгоритм Стефеста)
    td - безразмерное давление, число или numpy массив
    rd - безразмерный радиус, по умолчанию rd=1 - соответствует давлению на забое
         должно быть числом
    результат массив массивов давления от времени
    """
    pd_inv = anaflow.get_lap_inv(pd_lapl_line_source, rd=rd)
    return pd_inv(td)

# для сравнения приведем функции из лекции 1
# Решение линейного стока уравнения фильтрации с использованием Ei
def pd_ei(td, rd=1.):
    """
    Решение линейного стока уравнения фильтрации
    td - безразмерное время, число или numpy массив, больше нуля
    rd - безразмерное расстояние, по умолчанию rd=1 - соответствует давлению на забое
         число или numpy массив
    результат массив массивов давления от времени
    """
    return -0.5 * spec.expi(-0.25 * rd * rd / td)

Пример построения графиков сравнения реализаций решений линейного стока в безразмерных координатах

# построим график безразмерного давления от расстояния
# двумя способами и их разницу также
td = 100
rd_arr = np.logspace(1, 2 ,200)
pd_ei_arr = pd_ei(td, rd_arr) 
pd_lap_inv_arr = np.squeeze([pd_line_source_inv(td, rd_i) for rd_i in rd_arr])

# при построении используем векторный расчет
plt.rcParams["figure.figsize"] = (8,4)
fig, (ax1, ax2) = plt.subplots(1,2)

ax1.plot(rd_arr, pd_ei_arr)
ax1.plot(rd_arr, pd_lap_inv_arr)
ax1.set_title("Решение для td = {}".format(td))
ax1.set_xlabel("rd")
ax1.set_ylabel("pd")

ax2.plot(rd_arr, pd_ei_arr - pd_lap_inv_arr)
ax2.set_title("Разность решений")
ax2.set_xlabel("rd")
ax2.set_ylabel("pd")

plt.show()

Построенные графики показывают, что алгоритм Стефеста для численного обратного преобразования Лапласа дает погрешность до 2e-5 или примерно в пятом знаке после запятой, что является удовлетворительным результатом для большинства практических расчетов.

Определение функции для перевода размерных и безразмерных величин

# определим функции для перевода размерных переменных в безразмерные и обратно
# пригодится потом для построения графиков и ведения расчетов

# при наименовании функций придерживаемся следующих соглашений
# сначала идет название того, что считаем
# в конце указывается размерность результата, если это уместно

def r_from_rd_m(rd, rw_m=0.1):
    """
    перевод безразмерного расстояния в размерное
    rd -  безразмерное расстояние
    rw_m -  радиус скважины, м
    """
    return rd*rw_m

def rd_from_r(r_m, rw_m=0.1):
    """
    перевод размерного расстояния в безразмерное
    r_m - размерное расстояние, м
    rw_m - радиус скважины, м
    """
    return r_m/rw_m

def t_from_td_hr(td, k_mD=10, phi=0.2, mu_cP=1, ct_1atm=1e-5, rw_m=0.1):
    """
    перевод безразмерного времени в размерное, результат в часах
    td - безразмерное время
    k_mD - проницаемость пласта, мД
    phi - пористость, доли единиц
    mu_cP - динамическая вязкость флюида, сП
    ct_1atm - общая сжимаемость, 1/атм
    rw_m - радиус скважины, м
    """
    return td * phi * mu_cP * ct_1atm * rw_m * rw_m / k_mD / 0.00036

def td_from_t(t_hr, k_mD=10, phi=0.2, mu_cP=1, ct_1atm=1e-5, rw_m=0.1):
    """
    перевод размерного времени в безразмерное
    t_hr - размерное время, час
    k_mD - проницаемость пласта, мД
    phi - пористость, доли единиц
    mu_cP - динамическая вязкость флюида, сП
    ct_1atm - общая сжимаемость, 1/атм
    rw_m - радиус скважины, м
    """
    return  0.00036 * t_hr * k_mD / (phi * mu_cP * ct_1atm * rw_m * rw_m) 

def p_from_pd_atma(pd, k_mD=10, h_m=10, q_sm3day=20, b_m3m3=1.2, mu_cP=1, pi_atma=250):
    """
    перевод безразмерного давления в размерное, результат в абсолютных атмосферах
    pd - безразмерное давление
    k_mD - проницаемость пласта, мД
    h_m - мощность пласта, м
    q_sm3day - дебит на поверхности, м3/сут в с.у.
    fvf_m3m3 - объемный коэффициент нефти, м3/м3
    mu_cP - динамическая вязкость флюида, сП
    pi_atma - начальное давление, абсолютные атм.
    """
    return pi_atma - pd * 18.42 * q_sm3day * b_m3m3 * mu_cP / k_mD / h_m 

def pd_from_p(p_atma, k_mD=10, h_m=10, q_sm3day=20, b_m3m3=1.2, mu_cP=1, pi_atma=250):
    """
    перевод размерного давления в безразмерное
    p_atma - давление
    k_mD - проницаемость пласта, мД
    h_m - мощность пласта, м
    q_sm3day - дебит на поверхности, м3/сут в с.у.
    fvf_m3m3 - объемный коэффициент нефти, м3/м3
    mu_cP - динамическая вязкость флюида, сП
    pi_atma - начальное давление, абсолютные атм.
    """
    return (pi_atma - p_atma) / (18.42 * q_sm3day * b_m3m3 * mu_cP) * k_mD * h_m 

Пример построения графиков сравнения реализаций решений линейного стока в размерных координатах

# построим график безразмерного давления от расстояния в безразмерных переменных
# двумя способами и их разницу также
td = 100000
rd_arr = np.arange(1, 1000 ,10)
# при построении используем векторный расчет
pd_ei_arr = pd_ei(td, rd_arr) 
pd_lap_inv_arr = np.squeeze([pd_line_source_inv(td, rd_i) for rd_i in rd_arr])

plt.rcParams["figure.figsize"] = (8,4)
fig, (ax1, ax2) = plt.subplots(1,2)
ax1.plot(r_from_rd_m(rd_arr), p_from_pd_atma(pd_ei_arr) )
ax1.plot(r_from_rd_m(rd_arr), p_from_pd_atma(pd_lap_inv_arr) )

ax2.plot(r_from_rd_m(rd_arr), p_from_pd_atma(pd_ei_arr) - 
                              p_from_pd_atma(pd_lap_inv_arr))
ax1.set_title("Решение уравнения для t = {:.2f}".format(t_from_td_hr(td)))
ax2.set_title("Разность расчетов двумя способами")
ax1.set_xlabel("r, м")
ax1.set_ylabel("p, атм")
ax2.set_xlabel("r, м")
ax2.set_ylabel(r"$\Delta p$, атм")
plt.tight_layout()
plt.show()

При переводе значений давления в размерные величины можно получить, что ошибка в пятом знаке после запятой соответствует примерно 10 Па.

1.4 Решение для конечного радиуса скважины

Решение для конечного радиуса скважины в пространстве Лапласа будет иметь вид

\[ \tilde{p}_D(u) = \frac{K_0 \left( r_D \sqrt u \right)}{u \sqrt u K_1(\sqrt u)} \]

где

\(K_0\), \(K_1\) - модифицированные функции Бесселя

# пример функции реализующий расчет решения в пространстве Лапласа
def pd_lapl_finite_rw(u, rd=1.):
    """
    расчет решения c конечным радиусом скважины для безразмерного давления в пространстве Лапласа
    u - переменная пространства Лапласа
    rd- безразмерное расстояние от центра скважины
    """
    # полезно учесть, что при u>5e5 выражение kn(1, u05) обратится в ноль и будет деление на ноль
    # но если принудительно сделать там выражение равное нулу, обратное преобразование Лапласа
    # может выдавать очень странные результаты, поэтому лучше пока оставить как есть
    u05 = u**0.5
    return np.divide(spec.kn(0, rd * u05) , (u * u05 *  spec.kn(1, u05)))


# функция расчета безразмерного давления с использованием алгоритма Стефеста
# для численного обратного преобразования Лапласа
def pd_finite_rw_inv(td, rd=1.):
    """
    расчет решения c конечным радиусом скважины для безразмерного давления
    на основе численного обратного преобразования Лапласа (алгоритм Стефеста)
    td - безразмерное давление, число или numpy массив
    rd - безразмерный радиус, по умолчанию rd=1 - соответствует давлению на забое
    результат массив массивов давления от времени
    """
    pd_inv = anaflow.get_lap_inv(pd_lapl_finite_rw, rd=rd)
    return pd_inv(td)

1.5 Сравнение решений линейного стока и конечного радиуса скважины

Построим графики сравнения решений линейного стока и решения с учетом конечного радиуса скважины. Графики будем строить для фиксированного расстояния в зависимости от времени.

# построим график безразмерного давления от расстояния в безразмерных переменных
# двумя способами и их разницу также
td_arr = np.logspace(-4, 5 ,100)
rd = 1  # забойное давление

pd_finite_rw_inv_arr = pd_finite_rw_inv(td_arr, rd)
pd_line_source_inv_arr = pd_line_source_inv(td_arr, rd) 

# при построении используем векторный расчет
plt.rcParams["figure.figsize"] = (8,4)
fig, (ax1, ax2) = plt.subplots(1,2)
ax1.plot(t_from_td_hr(td_arr), p_from_pd_atma(pd_finite_rw_inv_arr) )
ax1.plot(t_from_td_hr(td_arr), p_from_pd_atma(pd_line_source_inv_arr) )
ax1.set_title("Сравнение решений")
ax1.set_xscale('log')
ax1.set_xlabel("t, час")
ax1.set_ylabel("p забойное, атм")

ax2.plot(t_from_td_hr(td_arr), p_from_pd_atma(pd_line_source_inv_arr) - 
                               p_from_pd_atma(pd_finite_rw_inv_arr))
ax2.set_title("Разность расчетов")
ax2.set_xscale('log')
ax2.set_xlabel("t, час")
ax2.set_ylabel(r"$\Delta p$, атм")
ax2.set_yscale('log')
plt.tight_layout()
plt.show()

Видно, что для забойного давления заметные расхождения наблюдаются для времен порядка 0.001 часа или 3 секунд. Следует помнить, что данная размерная оценка получена для заданных по умолчанию параметров скважин и изменение проницаемости или сжимаемости оценку изменит. Однако проведя эксперименты можно убедиться что время для которого расхождения будет заметно все равно останется относительно малым. В этом смысле оценка сходимости решений в безразмерных переменных предпочтительна (не будет зависеть от размерных парамтеров)

Посмотрим, что будет происходить вдали от скважины - меняем только параметр rd

# построим график безразмерного давления от расстояния в безразмерных переменных
# двумя способами и их разницу также
td_arr = np.logspace(1, 8 ,100)
rd = 1000  # забойное давление
# при построении используем векторный расчет
plt.rcParams["figure.figsize"] = (8,4)
fig, (ax1, ax2) = plt.subplots(1,2)
ax1.plot(t_from_td_hr(td_arr), p_from_pd_atma(pd_finite_rw_inv(td_arr, rd)) )
ax1.plot(t_from_td_hr(td_arr), p_from_pd_atma(pd_line_source_inv(td_arr, rd)) )
ax1.set_title("Сравнение решений на расстоянии {:.0f} метров".format(r_from_rd_m(rd)))
ax1.set_xscale('log')
ax1.set_xlabel("t, час")
ax1.set_ylabel("p, атм")

ax2.plot(t_from_td_hr(td_arr),p_from_pd_atma(pd_line_source_inv(td_arr, rd)) - 
                              p_from_pd_atma(pd_finite_rw_inv(td_arr, rd)) ) 
                               
ax2.set_title("Разность расчетов двумя способами")
ax2.set_xscale('log')
ax2.set_xlabel("t, час")
ax2.set_ylabel(r"$\Delta p$, атм")
plt.tight_layout()
plt.show()

Видно что вдали от скважины заметных расхождений в решениях нет. Поэтому для решения задач по анализу взаимного влияния скважин широкое применение получило относительно простое решение линейного стока.

2 Использование графиков с различными масштабами координат для анализа решений

При анализе гидродинамических исследований широко используются графики в полулогарифмических и двойных логарифмических координатах. Также часто возникает необходимость расчета производных и отображение их на графиках. Разберем технические вопросы построения таких графиков с использованием matplotlib.

Для примера будем использовать решения построенные ранее.

2.1 График в обычных координатах

Нас будут интересовать как графики от расстояния, так и графики от времени. Мы можем нарисовать несколько графиков одновременно с использованием matplotlib

# задаем параметры для отрисовки отдельных графиков
td = 10000
rd = 1
# задаем массивы для отрисовки графиков - количество точек на графике
rd_arr = np.logspace(1, 3, 100)
td_arr = np.logspace(1, 16, 100)

# создаем фигуру с двумя графиками, на котором все будет отрисовываться
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(8,5), gridspec_kw={'height_ratios': [3, 1]})

# задаем первый график используя оси ax1
ax1.plot(rd_arr, pd_ei(td, rd_arr) )
ax1.plot(-rd_arr, pd_ei(td, rd_arr) )
ax1.set(title="Воронка депресси для $t_D= {}$".format(td))
ax1.set_xlabel("$r_D$")
ax1.set_ylabel("$p_D$")

# задаем первый график используя оси ax2
ax2.plot(td_arr, pd_ei(td_arr, rd), color = 'red')
ax2.set(title="Динамика забойного давления")
ax2.set_xlabel("$t_D$")
ax2.set_ylabel("$p_D$")

ax3.axis('off')
ax4.plot([0,0,td_arr[-1]],[0, 1, 1], color = 'red')
ax4.set(title="Динамика дебита")
ax4.set_xlabel("$t_D$")
ax4.set_ylabel("$q_D$")
ax4.set_ylim(0, 1.5)

fig.tight_layout()
plt.show()

На построенных графиках мы видим решение в обычных - декартовых координатах. Такие данные мы получим при первичном анализе замерных данных с датчиков давления при запуске скважины (красный график). Для “ручного” анализа удобно преобразовать координаты использовав факт, что построенные решение похожи на логарифм.

2.2 График в полулогарифмических координатах для безразмерного решения

Графики похожие на логарифм при отображении в полулогарифмических координатах превращаются в прямую линию. Поэтому эти графики пользуются популярностью при анализе гидродинамических исследований.

Полулогарифмические координаты легко настраиваются всеми распространенными пакетами отображения. Ниже пример для matplotlib для решения с конечным радиусом скважины (12.22) построенным с использованием численного обратного преобразования Лапласа.

# задаем параметры для отрисовки отдельных графиков
td = 10000
rd = 1
# задаем массивы для отрисовки графиков - количество точек на графике
rd_arr = np.logspace(0, 3, 100)
td_arr = np.logspace(-3, 16, 100)

# создаем фигуру с двумя графиками, на котором все будет отрисовываться
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(8,5), gridspec_kw={'height_ratios': [3, 1]})

# задаем первый график используя оси ax1
ax1.plot(rd_arr, np.squeeze([pd_finite_rw_inv(td, rdi) for rdi in rd_arr]) )
ax1.plot(-rd_arr, np.squeeze([pd_finite_rw_inv(td, rdi) for rdi in rd_arr]) )
ax1.set(title="Воронка депресси для $t_D= {}$".format(td))
ax1.set_xlabel("$r_D$")
ax1.set_ylabel("$p_D$")
ax1.set_xscale('symlog', linthresh=1, linscale=0.6)

# задаем первый график используя оси ax2
ax2.plot(td_arr, pd_finite_rw_inv(td_arr, rd), color = 'red')
ax2.set(title="Динамика забойного давления")
ax2.set_xlabel("$t_D$")
ax2.set_ylabel("$p_D$")
ax2.set_xscale('log')

ax3.axis('off')
ax4.plot([0,0,td_arr[-1]],[0, 1, 1], color = 'red')
ax4.set(title="Динамика дебита")
ax4.set_xlabel("$t_D$")
ax4.set_ylabel("$q_D$")
ax4.set_ylim(0, 1.5)
ax4.set_xlim(ax2.get_xlim())
ax4.set_xscale('log')

fig.tight_layout()
plt.show()

Видно для небольших (относительно) расстояний и для больших (относительно) времен построенные решения выглядят как прямые линии, что означает, что они могут быть аппроксимированы с использованием логарифмического приближения.

Например решение линейного стока \[ p_D(r_D,t_D) = - \frac{1}{2} Ei \left(- \dfrac{ r_D^2}{4t_d} \right) +S \]

может быть описано следующим выражением

\[ p_D(r_D,t_D) = - \frac{1}{2} \ln \left( \dfrac{ r_D^2}{4t_d} \right) - \frac{1}{2}\gamma +S \]

\[ p_D(r_D,t_D) = \frac{1}{2}\ln(t_d) - \left[ \frac{1}{2} \ln \left( \dfrac{ r_D^2}{4} \right) + \frac{1}{2}\gamma -S\right] \] где \(\gamma = 0.57721566481\) - константа Эйлера

на графике от времени в полулогарифмических координатах логарифмическое приближение выглядит как кривая с наклоном \(0.5\)

2.3 График в полулогарифмических координатах для размерного решения

# задаем параметры для отрисовки отдельных графиков
t_hr = 10
# зададим два значения проницаемости что построить разные графики для сравнения
k_mD_1 = 10
k_mD_2 = 20
# определим дебит
q_scm3day = 25
# задаем массивы для отрисовки графиков - количество точек на графике
r_arr = np.logspace(-1, 3, 100)
t_arr = np.logspace(-7, 3, 100)

# создаем фигуру с двумя графиками, на котором все будет отрисовываться
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, 
                                             figsize=(8,5), 
                                             gridspec_kw={'height_ratios': [3, 1]})

# задаем первый график используя оси ax1
p_arr = lambda k: p_from_pd_atma(np.squeeze([pd_finite_rw_inv(td_from_t(t_hr), rdi) for rdi in rd_from_r(r_arr)]), k_mD=k, q_sm3day=q_scm3day)
ax1.plot(r_arr, p_arr(k_mD_1), label='k_mD = {}'.format(k_mD_1))
ax1.plot(r_arr, p_arr(k_mD_2), label='k_mD = {}'.format(k_mD_2))
ax1.plot(-r_arr, p_arr(k_mD_1), label='k_mD = {}'.format(k_mD_1))
ax1.plot(-r_arr, p_arr(k_mD_2), label='k_mD = {}'.format(k_mD_2))
ax1.set(title=f"Воронка депресси для $t= {t_hr:.2f}$ час")
ax1.set_xlabel("$r$")
ax1.set_ylabel("$p$")
ax1.set_xscale('symlog', linthresh=0.1, linscale=0.6)
ax1.legend()

# задаем первый график используя оси ax2
ax2.plot(t_arr, 
         p_from_pd_atma(pd_finite_rw_inv(td_from_t(t_arr), rd), k_mD=k_mD_1, q_sm3day=q_scm3day), color = 'red',
         label='k_mD = {}'.format(k_mD_1))
ax2.plot(t_arr, 
         p_from_pd_atma(pd_finite_rw_inv(td_from_t(t_arr), rd), k_mD=k_mD_2, q_sm3day=q_scm3day), color = 'magenta',
         label='k_mD = {}'.format(k_mD_2))
ax2.set(title="Динамика забойного давления")
ax2.set_xlabel("$t$")
ax2.set_ylabel("$p$")
ax2.set_xscale('log')
ax2.legend()

ax3.axis('off')
ax4.plot([0,0,t_arr[-1]],[0, q_scm3day, q_scm3day], color = 'red')
ax4.set(title="Динамика дебита")
ax4.set_xlabel("$t$")
ax4.set_ylabel("$q$")
ax4.set_ylim(0, q_scm3day*1.5)
ax4.set_xlim(ax2.get_xlim())
ax4.set_xscale('log')

fig.tight_layout()
plt.show()

На графиках в размерных переменных также радиальный приток (при котором справедливо логарифмическое приближение) отображается в виде прямых линий.

Решение в размерных переменных можно записать как \[ p\left(r,t\right)=p_i-\frac{18.42q_sB\mu}{kh}\left(-\frac{1}{2} Ei \left(-\frac{\varphi\mu c_tr^2}{0.00144kt}\right)+S\right) \]

Приведем несколько простых арифметических преобразований (чтобы не запутаться с коэффициентами) \[ p\left(r,t\right)=p_i-\frac{18.42q_sB\mu}{kh}\left(-\frac{1}{2} \ln \left(\frac{\varphi\mu c_tr^2}{0.00144kt}\right) -\frac{1}{2}\gamma + S\right) \]

\[ p\left(r,t\right)=p_i-\frac{18.42q_sB\mu}{kh}\left(\frac{1}{2} \ln \left(\frac{0.00144kt}{\varphi\mu c_tr^2}\right) -\frac{1}{2}\gamma + S\right) \]

\[ p\left(r,t\right)=p_i-\frac{9.205q_sB\mu}{kh}\left( \ln \left(\frac{kt}{\varphi\mu c_tr^2}\right) +ln(0.00144) -\gamma + 2S\right) \]

Наконец получим логарифмическое приближение \[ p\left(r,t\right)=p_i-\frac{9.205q_sB\mu}{kh}\left(ln{\frac{kt}{\varphi\mu c_tr^2}}-7.12+2S\right) \]

Или выделив \(\ln(t)\) в явном виде получим \[ p\left(r,t\right)=-\frac{9.205q_sB\mu}{kh}\ln(t) + \left[ p_i-\frac{9.205q_sB\mu}{kh}\left(ln{\frac{k}{\varphi\mu c_tr^2}}-7.12+2S\right) \right] \]

Из этого выражения видно, что определив наклон графика в полулогарифмических координатах и его смещение по вертикали мы сможем определить величины \(\dfrac{9.205q_sB\mu}{kh}\) и \(\left[ p_i-\dfrac{9.205q_sB\mu}{kh}\left(ln{\dfrac{k}{\varphi\mu c_tr^2}}-7.12+2S\right) \right]\) откуда чаще все находят \(kh\) и \(S\).

2.4 График в двойных логарифмических координатах

Графики в двойных логарифмических координатах также весьма популярны при интерпретации гидродинамической модели.

Такие графики обладают несколькими полезными свойствами. 1. Степенные зависимости будут переводиться в прямые линии. Функция вида \(y=x^a\) при логарифмировании обоих частей равенства преобразуется к виду \(ln(y) = a \cdot ln(x)\), что соответствует прямой линии на лог-лог графике с наклом \(a\). 2. Зависимости вида \(y = a \cdot f(b \cdot x)\) преобразуется при логарифмировании к виду $ln(y) = ln(a) + ln(f(e^{ln(b) + ln(x)})) $ или $ln(y) = ln(a) + g({ln(b) + ln(x)}) $ где \(g(z) = ln(f(e^z))\). Таким образом для произвольной функции при отображении ее в лог-лог координатах множители аргумента и самой функции превращаются в сдвиги некоторой функции \(g\), которая остается неизменной. За счет этого сравнивая разные изображения функций отличающихся множителями - их можно вычислить определяя величины сдвигов.

Ниже приведен пример отображения графиков степенных функций в обычных и двойных логарифмических координатах

# задаем массивы для отрисовки графиков - количество точек на графике
rr = np.arange(1, 10, 1)


# создаем фигуру с двумя графиками, на котором все будет отрисовываться
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(8,4))

# задаем первый график используя оси ax1
ax1.plot(rr, rr**2 , label = '$x^2$')
ax1.plot(rr, rr**0.5, label = '$x^{0.5}$' )
ax1.plot(rr, rr**1, label = '$x$' )
ax1.set(title="Степенные функции")
ax1.set_xlabel("$r_D$")
ax1.set_ylabel("$p_D$")
ax1.legend()


# задаем первый график используя оси ax2
ax2.plot(rr, rr**2, label = '$x^2$')
ax2.plot(rr, rr**0.5, label = '$x^{0.5}$' )
ax2.plot(rr, rr**1, label = '$x$' )
ax2.set(title="Степенные функции в лог-лог")
ax2.set_xlabel("$t_D$")
ax2.set_ylabel("$p_D$")
ax2.set_xscale('log')
ax2.set_yscale('log')
ax2.legend()

fig.tight_layout()
plt.show()

2.5 Расчет и отображение значений производной функции заданной по точкам

Совместно с двойными логарифмическими координатами часто применяется логарифмическая производная давления. Численный расчет производной несложен - например можно посчитать методом конечных разностей. При использовании numpy массивов можно воспользоваться функцией numpy.diff, которая вычисляет разницу между соседними элементами массива. Тут надо учитывать, что длина массива после дифференцирования будет меньше при использовании функции numpy.diff. Для того чтобы посчитать производную функции f(x) выражение с использованием np.diff будет иметь вид np.diff(f(x))/np.diff(x).

Другой вариант можно использовать функцию numpy.gradient куда надо передать табличные значения для дифференцирования и соответствующие значения аргумента.

При дифференцировании реальных данных полезно помнить, что операция численного дифференцирования может оказаться неустойчивой. В таких случаях необходимо применять разные виды фильтрации данных и стараться использовать более устойчивые алгоритмы численного дифференцирования (центральные разности с большим шагом и т.п.)

td = np.arange(1, 100, 0.1)
rr = 1

# при расчете разностей надо учесть, что длина массива уменьшится и убрать один элемент массива аргументов
plt.plot(td[:-1], np.diff(pd_ei(td))/np.diff(td), label='diff')
plt.plot(td, np.gradient(pd_ei(td), td), label='gradient')
plt.plot(td, pd_ei(td) , label='f')

plt.title("Пример функции и ее производной")
plt.xlabel("$t_D$")
plt.ylabel("$p_D$")
plt.legend()
plt.show()

Интересно отметить, что расчет numpy.gradient более универсален и проще в применении, однако простой расчет через разности немного быстрее.

%timeit np.diff(pd_ei(rr,td))/np.diff(td)
%timeit np.gradient(pd_ei(rr,td), td)

175 μs ± 2.2 μs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
213 μs ± 5.45 μs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

Приведем сравнение двух вариантов расчета производной - разностями и с использованием numpy.gradient для двух наборов данных - с небольшим количеством точек (10 точек) и большим (300) для логарифмически распределенных точек

td = np.logspace(1, 2, 10)
td_fine = np.logspace(1, 2, 300)
rr = 1

# при расчете разностей надо учесть, что длина массива уменьшится и убрать один элемент массива аргументов
plt.plot(td[:-1], np.diff(pd_ei(td))/np.diff(td), label='diff')
plt.plot(td, np.gradient(pd_ei(td), td, edge_order=2), label='gradient')
plt.plot(td_fine[:-1], np.diff(pd_ei(td_fine))/np.diff(td_fine), label='diff_fine')
plt.plot(td_fine, np.gradient(pd_ei(td_fine), td_fine, edge_order=2), label='gradient_fine')
#plt.plot(td, pd_ei(td) , label='f')

plt.title("Сравнение расчета производной разными функциями")
plt.xlabel("$t_D$")
plt.ylabel("$p_D$")
plt.xscale('log')
plt.legend()
plt.show()

Видно что расчет разностями для небольшого числа точек может давать большую погрешность в отличии от numpy.gradient, особенно при применении опции edge_order=2 для numpy.gradient позволяющей более точно оценить значение производной на границе.

Другой пример можно привести для неравномерное заданного диапазона аргументов (что удобно делать например в Excel). В этом случае расчет производной с использованием numpy.gradient также работает заметно лучше

td1 = list(range(1,10,1)) + list(range(10,100,10)) + list(range(100,1000,100))

td = np.array(td1)
td_fine = np.logspace(0, 2, 300)
rr = 1

# при расчете разностей надо учесть, что длина массива уменьшится и убрать один элемент массива аргументов
plt.plot(td[:-1], np.diff(pd_ei(td))/np.diff(td), label='diff uneven')
plt.plot(td, np.gradient(pd_ei(td), td, edge_order=2), label='gradient uneven')
plt.plot(td_fine[:-1], np.diff(pd_ei(td_fine))/np.diff(td_fine), label='diff_fine')
plt.plot(td_fine, np.gradient(pd_ei(td_fine), td_fine, edge_order=2), label='gradient_fine')

plt.title("Пример расчета производной для неравномерно заданного диапазона аргументов")
plt.xlabel("$t_D$")
plt.ylabel("$p_D$")
plt.xscale('log')
plt.yscale('log')
plt.legend()
plt.show()

Построим разные варианты графиков с отображением безразмерного давления и производной

td = np.logspace(0, 3, 100)
rr = 1
pd = pd_ei(td)
dpd_dtd = np.gradient(pd_ei(td), td, edge_order=2)
dpd_dlntd = np.gradient(pd_ei(td), np.log(td), edge_order=2) 

fig, ax = plt.subplots(3,2, figsize=(8,10))
ax[0,0].set_title("Изменение функции и ее производной")
ax[0,0].plot(td, pd_ei(td) )
ax[0,0].plot(td, dpd_dtd)
ax[0,0].set_xscale('log')
ax[0,0].set_xlabel("$t_D$")
ax[0,0].set_ylabel("$p_D$")

ax[0,1].set_title("Изменение функции и ее производной")
ax[0,1].plot(td, pd_ei(td) )
ax[0,1].plot(td, dpd_dtd)
ax[0,1].set_xscale('log')
ax[0,1].set_yscale('log')
ax[0,1].set_xlabel("$t_D$")
ax[0,1].set_ylabel("$p_D$")

ax[1,0].set_title("Изменение функции и ее логарифмической производной")
ax[1,0].plot(td, pd_ei(td) )
ax[1,0].plot(td, dpd_dlntd)
ax[1,0].set_xscale('log')
ax[1,0].set_xlabel("$t_D$")
ax[1,0].set_ylabel("$p_D$")

ax[1,1].set_title("Изменение функции и ее логарифмической производной")
ax[1,1].plot(td, pd_ei(td) )
ax[1,1].plot(td, dpd_dlntd)
ax[1,1].set_xscale('log')
ax[1,1].set_yscale('log')
ax[1,1].set_xlabel("$t_D$")
ax[1,1].set_ylabel("$p_D$")

ax[2,0].set_title("Изменение функции и ее производной")
ax[2,0].plot(td, pd_ei(td) )
ax[2,0].plot(td, dpd_dtd)
ax[2,0].set_xlabel("$t_D$")
ax[2,0].set_ylabel("$p_D$")

ax[2,1].set_title("Изменение функции и ее производной")
ax[2,1].plot(td, pd_ei(td) )
ax[2,1].plot(td, dpd_dlntd)
ax[2,1].set_xlabel("$t_D$")
ax[2,1].set_ylabel("$p_D$")

plt.tight_layout()
plt.show()

3 Учет скин-фактора и послепритока для нестационарного решения

3.1 Послеприток к скважине

Чаще всего управление дебитом скважины происходит на поверхности - на скважинной арматуре. Регулируя задвижки можно открыть - запустить скважину или остановить работу скважины. Замер дебита также проводится измерительным устройством на поверхности. При этом строя решения уравнения фильтрации мы описываем работу пласта и предполагаем, что знаем дебит именно пласта - дебит на границе соединения скважины с пластом. Объем ствола скважины и сжимаемость флюида в скважине могут привести к тому, что показатели притока на забое и на поверхности будут отличаться, особенно на переходных режимах работы. Именно этот эффект учитывается моделями ствола скважины или моделя послепритока (wellbore storage or afterflow).

Самая простая модель предполагает постоянную сжимаемость ствола скважины. Для ее описания вводится коэффициент влияния ствола скважины

\[ C=V_wc_w = \dfrac{\Delta V}{\Delta P} \]

где

• \(C\) - коэффициент влияния ствола скважины, м\(^3\)/атм
• \(V_w\) - объем ствола скважины, м\(^3\)
• \(c_w\) - сжимаемость флюда (газожидкостной смеси) в стволе скважины, 1/атм
• \(\Delta V\) - изменение объема флюида в скважине, м\(^3\)
• \(\Delta P\) - изменение давления в скважина, атм

Коэффициент ствола скважины легко выразить для некоторых случаев.

• Для нагнетательной скважины

\[ C=V_wc_w \]

• Для фонтанирующей нефтяной скважины

\[ C=V_lc_l + \dfrac{ V_g}{P_g} \]

• Для механизированной скважины с динамическим уровнем

\[ = f\dfrac{A}{\rho g} \]

здесь

• \(C\) - коэффициент влияния ствола скважины, м\(^3\)/атм
• \(V_w\) - объем ствола скважины, м\(^3\)
• \(c_w\) - сжимаемость воды в стволе скважины, 1/атм
• \(V_l\) - объем жидкости в стволе фонтанирующей скважины, м\(^3\)
• \(c_l\) - сжимаемость жидкости в стволе фонтанирующей скважины, 1/атм
• \(V_g\) - объем газа в стволе фонтанирующей скважины, м\(^3\)
• \(P_g\) - давление газа в стволе фонтанирующей скважины, атм
• \(A\) - площадь внутреннего сечения ствола скважины, м\(^2\)
• \(\rho\) - плотность жидкости в стволе скважины, кг/м\(^3\)
• \(g\) - ускорение свободного падения, м/сек\(^2\)
• \(f\) - переводной коэффициент для согласования размерностей \(f=101325\)

Используя приведенные соотношения можно оценить величины коэффициента послепритока для разных случаев

• Для нагнетательной скважины

\[ C=V_wc_w = 30 \cdot 4.5\cdot10^{-5} = 1.35\cdot10^{-5} [м^3/атм] \]

• Для фонтанирующей нефтяной скважины

\[ C=V_lc_l + \dfrac{ V_g}{P_g} = 29 \cdot 5\cdot10^{-5} + 1 \cdot \dfrac{1}{100} = 0.01 [м^3/атм] \]

• Для механизированной скважины с динамическим уровнем

\[ C = f\dfrac{A}{\rho g} = \dfrac{101325 \cdot 0.015}{9.8 \cdot 800} = 0.2 [м^3/атм] \]

Для связи дебита скважины на поверхности \(q_{wh}\) и на забое \(q_{sf}\) можно получить выражение

\[ q_{sf} = q_{wh} B + 24 C \dfrac{dP_{wf}}{dt} \]

где

• \(q_{sf}\) - дебит на забое скважины, м\(^3\)/сут
• \(q_{wh}\) - дебит на устье скважины м\(^3\)/сут
• \(dP_{wf}\) - изменение забойного давления, атм
• \(dt\) - изменение времени, час
• \(C\) - коэффициент влияния ствола скважины, м\(^3\)/атм
• \(B\) - объемный коэффициент нефти

4 Построение решения в пространстве Лапласа для решения с учетом послепритока и скин-фактора

Рассматривается уравнение фильтрации в безразмерных переменных

\[ \frac{\partial p_D}{ \partial t_D} = \frac{1}{r_D}\left[ \frac{ \partial{}}{ \partial{r_D} }\left( r_D \dfrac{\partial p_D}{ \partial r_D} \right) \right] \]

вводятся следующие безразмерные переменные

• \(r_D\) - безразмерное расстояние от центра скважины
• \(t_D\) - безразмерное время
• \(p_D\) - безразмерное давление

Соответствующие определения безразмерных переменных

\[ r_D = \frac{r}{r_w} \] \[ t_D = \frac{0.00036 kt}{\phi \mu c_t r_w^2} \] \[ p_D = \frac{kh}{ 18.42 q B \mu} \left( p_i - p \right) \]

где в свою очередь

• \(q\) - дебит скважины на поверхности, приведенный к нормальным условиям, м3/сут
• \(\phi\) - пористость, доли единиц
• \(\mu\) - вязкость нефти в пласте, сП
• \(B\) - объемный коэффициент нефти, м3/м3
• \(p_i\) - начальное давление в пласте, атм
• \(p\) - давление забойное, атм
• \(c_t\) - общая сжимаемость системы в пласте, 1/атм
• \(k\) - проницаемость, мД
• \(t\) - время, час
• \(r\) - расстояние от центра скважины, м
• \(r_w\) - радиус скважины, м

дополнительно вводим

\[ q_D = \frac{ q_{sf} }{qB} \]

• \(q_D\) - безразмерный дебит на поверхности. Нормириуется на произвольный дебит \(q\)
• \(q_{sf}\) - дебит на поверхности. Может меняться для задачи постоянного дебита из пласта

\[ C_D = 0.159 \dfrac{C}{h\phi \mu c_t r_w^2} \]

• \(C_D\) - безразрмерный коэффициент влияния скважины

4.1 Граничные условия для конечного радиуса скважины

• начальное условие. До запуска скважины в момент времени \(t_D = 0\) давление в пласте равно начальному во всех точках \(p=p_i\) \[ t_D < 0, p_D = 0 \]

• условие постоянства дебита на скважине - граничное условие на скважине

\[ \lim_{r_D \to 1} {r_D \frac{\partial p_D}{\partial r_D}} = -1 \]

• условие на бесконечном расстоянии возмущения от скважине нет \[ r_D = \infty, p_D = 0 \]

Общее решение имеет вид

\[ \tilde{p}_D(u, r_D) = A(u) K_0(r_D \sqrt u) + B(u) I_0(r_D \sqrt u) \]

Частное решение соответсвующее приведенным условиям будет иметь вид

\[ \tilde{p}_D(s) = \frac{1}{s^{\frac{3}{2}}} \frac{ K_0 \left( r_D \sqrt s \right) }{K_1 \left( \sqrt s \right)} \]

где

\(K_0\), \(K_1\) - модифицированные функции Бесселя

4.2 Граничные условия для скина и послепритока

• начальное условие. До запуска скважины в момент времени \(t_D = 0\) давление в пласте равно начальному во всех точках \(p=p_i\) \[ t_D < 0, p_D = 0 \]

• условие на бесконечном расстоянии возмущения от скважине нет \[ r_D = \infty, p_D = 0 \]

• условие постоянства дебита на скважине - граничное условие на скважине \[ \lim_{r_D \to 1} {r_D \frac{\partial p_D}{\partial r_D}} = -q_D \]

где

\[ q_D = 1-C_D \dfrac{dp_{wfD}}{dt_D} \]

где \[ p_{wfD} = p_D + q_D S \]

Общее решение имеет вид

\[ \tilde{p}_D(u, r_D) = A(u) K_0(r_D \sqrt u) + B(u) I_0(r_D \sqrt u) \]

Частное решение решение с учетом скин-фактор и послепритока можно представить в виде

\[ \tilde{p}_{wbD}(s) = \frac{1}{s} \left[ \frac{S + s \tilde{p}_D(r_d=1,s) }{1 + s C_D (S + s \tilde{p}_D(r_d=1,s))} \right] \]

где \(\tilde{p}_D(r_d=1,s)\) - известное частное решение

4.3 Расчет на python и построение графиков

# решение линейного стока в пространстве Лапласа
def pd_lapl_wbs(pd_lapl, S=0, Cd=0):
    def pdl(u, rd, S=0, Cd=0):
        upd_lapl = u * pd_lapl(u, rd=rd)
        return 1 / u * (S + upd_lapl) / (1 + u * Cd * (S + upd_lapl))
    return pdl

# построим функцию инвертирующую решение в пространстве Лапласа

# функция расчета безразмерного давления с использованием алгоритма Стефеста
# для численного обратного преобразования Лапласа
def pd_inv(pd_lapl, td, rd, **kwargs):
     
    
    if isinstance(rd, np.ndarray):
        # если на входе массив, то вручную итерируем по всем элементам и 
        # собираем массив массивов
        return np.array(list(map(lambda r: pd_inv(pd_lapl, td, rd=r, **kwargs), rd)))
    else:
        pd_inv_ = anaflow.get_lap_inv(pd_lapl, rd=rd, **kwargs)
        return pd_inv_(td)

t_d = np.arange(1, 1e8, 1e6)
r_d = 1
S = 0
C_d = 100000 
p_d_1 = pd_inv(pd_lapl_finite_rw, td=t_d, rd=r_d)


fig, ax = plt.subplots(1,1)

clist = np.arange(1,10,1)
for ci in clist:
    p_d_2 = pd_inv(pd_lapl_wbs(pd_lapl_finite_rw), td=t_d, rd=r_d, S=S, Cd=ci*C_d)
    plt.plot(t_d, p_d_2)

#ax.set_xscale('log')
plt.plot(t_d, p_d_1)
plt.show()

t_d = np.logspace(1, 10)
r_d = 1
S = 0
C_d = 10000 
p_d_1 = pd_inv(pd_lapl_finite_rw, td=t_d, rd=r_d)

fig, ax = plt.subplots(1,1)

clist = np.arange(1,100,10)
for ci in clist:
    p_d_2 = pd_inv(pd_lapl_wbs(pd_lapl_finite_rw), td=t_d, rd=r_d, S=S, Cd=ci*C_d)
    plt.plot(t_d, p_d_2)

ax.set_xscale('log')
plt.plot(t_d, p_d_1)
plt.show()

t_d = np.logspace(1, 10)
r_d = 1
S = 0
C_d = 10000 
p_d_1 = pd_inv(pd_lapl_finite_rw, td=t_d, rd=r_d)


fig, ax = plt.subplots(1,1)

clist = np.arange(1,100,10)
for ci in clist:
    p_d_2 =  pd_inv(pd_lapl_wbs(pd_lapl_finite_rw), td=t_d, rd=r_d, S=S, Cd=ci*C_d)
    plt.plot(t_d, p_d_2)

ax.set_xscale('log')
ax.set_yscale('log')
plt.plot(t_d, p_d_1)
plt.show()

5 Уравнение фильтрации и его решения - построение с использованием python и sympy

5.1 Уравнение фильтрации

Рассматривается уравнение фильтрации в безразмерных переменных

\[ \frac{\partial p_D}{ \partial t_D} = \frac{1}{r_D}\left[ \frac{ \partial{}}{ \partial{r_D} }\left( r_D \dfrac{\partial p_D}{ \partial r_D} \right) \right] \]

где введены следующие безразмерные переменные * \(r_D\) - безразмерное расстояние от центра скважины * \(t_D\) - безразмерное время * \(p_D\) - безразмерное давление

Соответствующие определения безразмерных переменных

\[ r_D = \frac{r}{r_w} \] \[ t_D = \frac{0.00036 kt}{\phi \mu c_t r_w^2} \]

\[ p_D = \frac{kh}{ 18.42 q_s B \mu} \left( p_i - p \right) \]

где в свою очередь

• \(q_s\) - дебит скважины на поверхности, приведенный к нормальным условиям, м3/сут
• \(\phi\) - пористость, доли единиц
• \(\mu\) - вязкость нефти в пласте, сП
• \(B\) - объемный коэффициент нефти, м3/м3
• \(p_i\) - начальное давление в пласте, атм
• \(p\) - давление забойное, атм
• \(c_t\) - общая сжимаемость системы в пласте, 1/атм
• \(k\) - проницаемость, мД
• \(t\) - время, час
• \(r\) - расстояние от центра скважины, м
• \(r_w\) - радиус скважины, м

Решение такого уравнение может быть получено с использованием преобразования Лапласа

\[ L \left [ f(t) \right] = \tilde{f}(u) = \int_{0}^{\infty}f(t)e^{-ut}dt \]

где \(u\) параметр пространства Лапласа соответствующий времени.

Тогда уравнение в пространстве Лапласа преобразуется к виду:

\[ u \tilde{p}_D = \dfrac{1}{r_D} \left[\dfrac{d}{d r_D} \left(r_D \dfrac{d{\tilde{p}_D}}{d r_D} \right) \right] \tag{5.1}\]

Уравнение (5.1) известно как модифицированное уравнение Бесселя где

• \(u\) - переменная пространства Лапласа, соответствующая времени
• \(\tilde{p}_D(u, r_D)\) - изображение давления в пространстве Лапласа

С уравнением (5.1) можно начать работать с использованием системы компьютерной алгебры sympy.

# определим переменные с которыми будем работать 

r_d = sp.symbols('r_d', real=True, positive=True)
p_d = sp.symbols('p_d', cls=sp.Function, real=True, positive=True)
u = sp.symbols('u',real=True, positive=True)
C1, C2 = sp.symbols('C1 C2')

# определим уравнение фильтрации в пространстве Лапласа
diff_eq = sp.Eq(u * p_d(r_d),
                1 / r_d * (sp.diff(r_d * sp.diff(p_d(r_d), r_d) , r_d)) )
diff_eq

\(\displaystyle u p_{d}{\left(r_{d} \right)} = \frac{r_{d} \frac{d^{2}}{d r_{d}^{2}} p_{d}{\left(r_{d} \right)} + \frac{d}{d r_{d}} p_{d}{\left(r_{d} \right)}}{r_{d}}\)

5.2 Общее решение модифицированного уравнения Бесселя

# попробуем решить уравнение 
soln = sp.dsolve(diff_eq, p_d(r_d))
soln

\(\displaystyle p_{d}{\left(r_{d} \right)} = C_{1} I_{0}\left(r_{d} \sqrt{u}\right) + C_{2} Y_{0}\left(i r_{d} \sqrt{u}\right)\)

# проверим, что решение полученное dsolve удовлетворяет уравнению
sp.checkodesol(diff_eq, soln)

(True, 0)

Уравнение Бесселя легко решается, \[ \displaystyle \operatorname{p_{d}}{\left(r_{d} \right)} = C_{1} I_{0}\left(r_{d} \sqrt{u}\right) + C_{2} Y_{0}\left(i r_{d} \sqrt{u}\right) \tag{5.2}\]

хотя решение (5.2) по виду и немного отличается от того, что обычно приводится в книгах. Вместо мофицированной функции Бесселя второго рода \(K_0(x)\) решение выражается через функцию Бесселя второго рода \(Y_0(ix)\) для мнимого арумента.

Известны выражения связывающие функции Бесселя первого \(J_\alpha(iz)\) и второго рода \(Y_\alpha(iz)\) для мнимых аргументов с модифицированными функциями Бесселя первого \(I_\alpha(z)\) и второго рода \(K_\alpha(z)\).

\[ J_\alpha(iz) = e^{\frac{\alpha\pi i}{2}} I_\alpha(z) \] \[ Y_\alpha(iz) = e^{\frac{(\alpha+1)\pi i}{2}}I_\alpha(z) - \frac{2}{\pi}e^{-\frac{\alpha\pi i}{2}}K_\alpha(z) \tag{5.3}\]

преобразуем (5.3)

\[ Y_0(iz) = e^{\frac{\pi i}{2}}I_0(z) - \dfrac{2}{\pi}K_0(z) \]

Учитывая эти выражения можно убедиться что выражение вида

\[ \tilde{p}_D(u, r_D) = A(u) K_0(r_D \sqrt u) + B(u) I_0(r_D \sqrt u) \tag{5.4}\]

которое можно найти в книгах также будет являться решение уравнения фильтрации, что можно проверить командой sympy.checkodesol

# зададим в явном виде решение с использованием K_0 и I_0
A, B = sp.symbols('A B')
soln2 = sp.Eq(p_d(r_d) , A * sp.besselk(0, r_d * sp.sqrt(u)) + B * sp.besseli(0, r_d * sp.sqrt(u)))
soln2

\(\displaystyle p_{d}{\left(r_{d} \right)} = A K_{0}\left(r_{d} \sqrt{u}\right) + B I_{0}\left(r_{d} \sqrt{u}\right)\)

# проверим, что это решение также удовлетворяет исходному уравнению
sp.checkodesol(diff_eq, soln2)

(True, 0)

Дальше покажем как можно работать с решением (5.4) для построения частных решений.

5.3 Частное решение для бесконечного пласта и конечного радиуса скважины.

Частное решение строится за счет поиска параметров \(A\) и \(B\) для уравнения (2.5) так, чтобы они удовлетворяли граничным условиям.

1. Граничное условия на бесконечности

\[ \lim_{r_D \to \infty} p_D(r_D, t_D) = 0 \]

в пространстве Лапласа преобразуется в следующее

\[ \lim_{r_D \to \infty} \tilde{p}_D(r_D, u) = 0 \]

2. Граничное условие на скважине

\[ \lim_{r_D \to r_{wD}} \left[ r_D \dfrac{ \partial p_D(r_D, t_D)}{\partial r_D} \right] = -1 \]

в пространстве Лапласа с учетом выражения (1.9) преобразуется к

\[ \lim_{r_D \to r_{wD}} \left[ r_D \dfrac{ d \tilde{p}_D}{d r_D} \right] = -\dfrac{1}{u} \]

Условие на бесконечности трудно записать в sympy в явном виде. Но можно проверить пределы выражения (2.5) при разных значениях параметров \(A\) и \(B\)

sp.limit(soln2.rhs, r_d, sp.oo)

\(\displaystyle \infty \operatorname{sign}{\left(B \right)}\)

Предел при призвольных \(A\) и \(B\) не определен.

Положим \(A=0\). Это делается командой .subs(A, 0) для исследуемого выражения

sp.limit(soln2.rhs.subs(A, 0), r_d, sp.oo)

\(\displaystyle \infty \operatorname{sign}{\left(B \right)}\)

Получаем, что этот предел также не определен. Положим \(B=0\).

sp.limit(soln2.rhs.subs(B, 0), r_d, sp.oo)

\(\displaystyle 0\)

Получаем, что при \(B=0\) выражение \(\displaystyle \lim_{r_{d} \to \infty}\left(A K_{0}\left(r_{d} \sqrt{u}\right) + B I_{0}\left(r_{d} \sqrt{u}\right)\right) = 0\)

Откуда можно сделать вывод, что для приведенных граничных условий \(B=0\).

soln3 = soln2.subs(B, 0)
soln3

\(\displaystyle p_{d}{\left(r_{d} \right)} = A K_{0}\left(r_{d} \sqrt{u}\right)\)

Найдем констранту \(A\) из второго граничного условия

\[ \lim_{r_D \to r_{wD}} \left[ r_D \dfrac{ d \tilde{p}_D}{d r_D} \right] = -\dfrac{1}{u} \]

bc2 = sp.Eq(r_d * soln3.rhs.diff(r_d), -1/u)
bc2

\(\displaystyle - A r_{d} \sqrt{u} K_{1}\left(r_{d} \sqrt{u}\right) = - \frac{1}{u}\)

r_wd = sp.symbols('r_wd', real=True, positive=True)
bc2_sol = sp.solve(bc2.subs(r_d, r_wd), A)
bc2_sol[0]

\(\displaystyle \frac{1}{r_{wd} u^{\frac{3}{2}} K_{1}\left(r_{wd} \sqrt{u}\right)}\)

soln4 = soln3.subs(A, bc2_sol[0])
soln4

\(\displaystyle p_{d}{\left(r_{d} \right)} = \frac{K_{0}\left(r_{d} \sqrt{u}\right)}{r_{wd} u^{\frac{3}{2}} K_{1}\left(r_{wd} \sqrt{u}\right)}\)

При \(r_{wd} = 1\) получим

soln5 = soln4.subs(r_wd, 1)
soln5

\(\displaystyle p_{d}{\left(r_{d} \right)} = \frac{K_{0}\left(r_{d} \sqrt{u}\right)}{u^{\frac{3}{2}} K_{1}\left(\sqrt{u}\right)}\)

Получили решение аналогичное решению (1.17). При этом часть преобразований при получении решений была сделана с использованием sympy.

\[ \displaystyle \operatorname{p_{d}}{\left(r_{d} \right)} = \frac{K_{0}\left(r_{d} \sqrt{u}\right)}{u^{\frac{3}{2}} K_{1}\left(\sqrt{u}\right)} \]

5.4 Построение численных решений на основе аналитического с использованием mpmath

Полученное решение можно использовать для отображения. Оно легко копируется в формате LaTeX и встравляется в блоки описания. Его можно использовать для дальнейших аналитических преобразований. Его можно также использовать для построения численных решений.

Покажем как решение можно построить с использованием mpmath и реализации численного обратного преобразования mpmath

# преобразование sympy выражения в lambda функцию с использованием mpmath
soln5_mpmath = sp.lambdify([u, r_d], soln5.rhs, modules = ['mpmath'])

# после преобразования lambda функцию можно использовать 
# для проведения расчетов как обычныую функцию

# пример расчета для одного значения, чтобы убедиться, что все работает
soln5_mpmath(u=1, r_d=1)

mpf('0.69948393559377231')

# определим функцию для расчета обратного преобразования Лапласа
# отдельная функция удобна чтобы передавать r_d второй параметр
def soln5_inv_lapl(t_d=1000, r_d=1):
    sol_lap = lambda u: soln5_mpmath(u=u, r_d=r_d)
    return mp.invertlaplace(sol_lap, t_d)

# пример расчета безразмерного давления с использованием mpmath
soln5_inv_lapl(t_d=1000, r_d=1)

mpf('3.8605905955786238')

# нарисуем график зависимости безразмерного давления от безразмерного времени
# поскольку считает относительно медленно ограничим количество точек для расчета
mp.plot(soln5_inv_lapl,[1,100], points=20)

# нарисуем график зависимости безразмерного давления от безразмерного расстояния
# поскольку считает относительно медленно ограничим количество точек для расчета
td=1000
mp.plot(lambda rd: soln5_inv_lapl(td,rd),[1,100], points=20)

5.5 Построение численных решений на основе аналитического с использованием numpy и scipy

Использованием mpmath удобно для проверки расчетов, так как эта библиотека позволяет выбрать разные методы расчета и задать нужную точность. Но расчеты при этом проводятся медленно. Для более быстрых расчетов можно использовать расчетные функции numpy и scipy, а для обратного преобразования Лапласа библиотеку anaflow

# преобразование sympy выражения в lambda функцию с использованием numpy
soln5_numpy = sp.lambdify([u, r_d], soln5.rhs, 
                           modules = ['numpy', 'scipy'])

# после преобразования lambda функцию можно использовать 
# для проведения расчетов как обычныую функцию

# пример расчета для одного значения, чтобы убедиться, что все работает
soln5_numpy(u=1, r_d=1)

np.float64(0.6994839355937723)

# определим функцию для расчета обратного преобразования Лапласа
# отдельная функция удобно чтобы передавать r_d второй параметр
def soln5_inv_lapl_numpy(t_d=1000, r_d=1):
    sol_lap = anaflow.get_lap_inv(soln5_numpy,  r_d=r_d)
    return sol_lap(t_d)

# пример расчета для одного значения, чтобы убедиться, что все работает
soln5_inv_lapl_numpy(1000)

array([3.86059079])

Построим графики с использованием сгенерированного решения. В отличии от mpmath решения считает быстрее, поэтому можно построить сразу много графиков.

# зададим диапазоны изменения параметров для построения графиков
td_arr = np.logspace(1,5,100)
rd_arr = np.logspace(1,3,100)

# построим графики
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[15,5])
for ri in np.linspace(1,200,10):
    ax1.plot(td_arr, soln5_inv_lapl_numpy(td_arr, r_d=ri), label = f'$r_d={ri:.0f}$')
for ti in np.linspace(1e4, 1e7,10):
    ax2.plot(rd_arr, np.squeeze( [soln5_inv_lapl_numpy(ti, r_d=ri) for ri in rd_arr ] ), label = f'$t_d={ti:.0f}$')
# наведем немного красоты
ax1.set_xlabel('$t_d$')
ax1.set_ylabel('$p_d$')
ax1.legend()
ax2.set_xlabel('$r_d$')
ax2.set_ylabel('$p_d$')
ax2.legend()
plt.show()

5.6 Частное решение для конечного круглого пласта с постоянным давлением на контуре и конечного радиуса скважины.

Частное решение строится за счет поиска параметров \(A\) и \(B\) для уравнения (5.4) так, чтобы они удовлетворяли граничным условиям.

1. Граничное условия на внешнем контуре с постоянным давлением

\[ p_D(r_e, t_D) = 0 \]

в пространстве Лапласа преобразуется в следующее

\[ \tilde{p}_D(r_e, u) = 0 \]

2. Граничное условие на скважине

\[ \lim_{r_D \to r_{wD}} \left[ r_D \dfrac{ \partial p_D(r_D, t_D)}{\partial r_D} \right] = -1 \]

в пространстве Лапласа преобразуется к

\[ \lim_{r_D \to r_{wD}} \left[ r_D \dfrac{ d \tilde{p}_D}{d r_D} \right] = -\dfrac{1}{u} \]

soln2

\(\displaystyle p_{d}{\left(r_{d} \right)} = A K_{0}\left(r_{d} \sqrt{u}\right) + B I_{0}\left(r_{d} \sqrt{u}\right)\)

r_e = sp.symbols('r_e', positive=True)
eq1 = sp.Eq(soln2.subs(r_d, r_e).rhs,0)
eq1

\(\displaystyle A K_{0}\left(r_{e} \sqrt{u}\right) + B I_{0}\left(r_{e} \sqrt{u}\right) = 0\)

eq2 = sp.Eq(r_d * soln2.rhs.diff(r_d), -1/u).subs(r_d, 1)
eq2

\(\displaystyle - A \sqrt{u} K_{1}\left(\sqrt{u}\right) + B \sqrt{u} I_{1}\left(\sqrt{u}\right) = - \frac{1}{u}\)

sol_AB = sp.solve([eq1,eq2], [A,B])
display(sol_AB[A])
display(sol_AB[B])

\(\displaystyle \frac{I_{0}\left(r_{e} \sqrt{u}\right)}{u^{\frac{3}{2}} I_{0}\left(r_{e} \sqrt{u}\right) K_{1}\left(\sqrt{u}\right) + u^{\frac{3}{2}} I_{1}\left(\sqrt{u}\right) K_{0}\left(r_{e} \sqrt{u}\right)}\)

\(\displaystyle - \frac{K_{0}\left(r_{e} \sqrt{u}\right)}{u^{\frac{3}{2}} I_{0}\left(r_{e} \sqrt{u}\right) K_{1}\left(\sqrt{u}\right) + u^{\frac{3}{2}} I_{1}\left(\sqrt{u}\right) K_{0}\left(r_{e} \sqrt{u}\right)}\)

sol_re_const_p = soln2.subs(A, sol_AB[A]).subs(B, sol_AB[B]).simplify()  
sol_re_const_p

\(\displaystyle p_{d}{\left(r_{d} \right)} = \frac{- I_{0}\left(r_{d} \sqrt{u}\right) K_{0}\left(r_{e} \sqrt{u}\right) + I_{0}\left(r_{e} \sqrt{u}\right) K_{0}\left(r_{d} \sqrt{u}\right)}{u^{\frac{3}{2}} \left(I_{0}\left(r_{e} \sqrt{u}\right) K_{1}\left(\sqrt{u}\right) + I_{1}\left(\sqrt{u}\right) K_{0}\left(r_{e} \sqrt{u}\right)\right)}\)

Таким образом получено решение для постоянного давления на круговой границе вокруг скважины \[ \displaystyle \operatorname{p_{d}}{\left(r_{d} \right)} = \frac{- I_{0}\left(r_{d} \sqrt{u}\right) K_{0}\left(r_{e} \sqrt{u}\right) + I_{0}\left(r_{e} \sqrt{u}\right) K_{0}\left(r_{d} \sqrt{u}\right)}{u^{\frac{3}{2}} \left(I_{0}\left(r_{e} \sqrt{u}\right) K_{1}\left(\sqrt{u}\right) + I_{1}\left(\sqrt{u}\right) K_{0}\left(r_{e} \sqrt{u}\right)\right)} \]

Ипользуя подход показанный выше можно построить графики зависимости безразмерного давления от времени и расстояния для этого решения с использованием numpy и scipy функций

# преобразование sympy выражения в lambda функцию с использованием numpy
soln_re_const_p_numpy = sp.lambdify([u, r_d, r_e], sol_re_const_p.rhs, 
                                     modules = ['numpy', 'scipy'])

# после преобразования lambda функцию можно использовать 
# для проведения расчетов как обычныую функцию

# определим функцию для расчета обратного преобразования Лапласа
# отдельная функция удобно чтобы передавать r_d второй параметр
def soln_re_const_p_inv_lapl_numpy(t_d=1000, r_d=1, r_e=100):
    sol_lap = anaflow.get_lap_inv(soln_re_const_p_numpy,  r_d=r_d, r_e=r_e)
    return sol_lap(t_d)

# зададим диапазоны изменения параметров для построения графиков
td_arr = np.logspace(1,5,100)
rd_arr = np.logspace(1,2,100)
re=100
# построим графики
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[15,5])
for ri in np.linspace(1,100,10):
    ax1.plot(td_arr, soln_re_const_p_inv_lapl_numpy(td_arr, r_d=ri, r_e=re), label = f'$r_d={ri:.0f}$')
for ti in np.linspace(1e2, 1e4,10):
    ax2.plot(rd_arr, np.squeeze( [soln_re_const_p_inv_lapl_numpy(ti, r_d=ri, r_e=re) for ri in rd_arr ] ), label = f'$t_d={ti:.0f}$')
# наведем немного красоты
ax1.set_xlabel('$t_d$')
ax1.set_ylabel('$p_d$')
ax1.legend()
ax2.set_xlabel('$r_d$')
ax2.set_ylabel('$p_d$')
ax2.legend()
plt.show()

5.6.1 Построение решения с использованием sympy.dsolve

Отметим, что в некоторых случаях можно построить частные решения для ОДУ напряму с использованием sympy

Вспомним с каким уравнением мы работаем

diff_eq

\(\displaystyle u p_{d}{\left(r_{d} \right)} = \frac{r_{d} \frac{d^{2}}{d r_{d}^{2}} p_{d}{\left(r_{d} \right)} + \frac{d}{d r_{d}} p_{d}{\left(r_{d} \right)}}{r_{d}}\)

Построим снова общее решение, чтобы убедиться, что на текущем этапе все работает

soln = sp.dsolve(diff_eq, p_d(r_d))
soln

\(\displaystyle p_{d}{\left(r_{d} \right)} = C_{1} I_{0}\left(r_{d} \sqrt{u}\right) + C_{2} Y_{0}\left(i r_{d} \sqrt{u}\right)\)

Построим частное решение сразу задав в sympy.dsolve граничные условия на искомую функцию

soln_sympy = sp.dsolve(diff_eq, p_d(r_d), 
                       ics={p_d(r_d).subs(r_d, r_e):0, 
                            p_d(r_d).diff(r_d).subs(r_d,1):-1/u}).simplify()
soln_sympy

\(\displaystyle p_{d}{\left(r_{d} \right)} = \frac{i \left(I_{0}\left(r_{d} \sqrt{u}\right) Y_{0}\left(i r_{e} \sqrt{u}\right) - I_{0}\left(r_{e} \sqrt{u}\right) Y_{0}\left(i r_{d} \sqrt{u}\right)\right)}{u^{\frac{3}{2}} \left(I_{0}\left(r_{e} \sqrt{u}\right) Y_{1}\left(i \sqrt{u}\right) - i I_{1}\left(\sqrt{u}\right) Y_{0}\left(i r_{e} \sqrt{u}\right)\right)}\)

Решение построено, хотя оно выражается через \(Y_0\) вместо \(K_0\) и в нем в явном виде присутствует \(i\) мнимая единица.

Это решение также можно использовать для проведения расчетов и построения графиков. Для того, чтобы подавить предупреждения об игнорировании мнимой части - при проведении обратного преобразования Лапласа оставим только реальную часть в явном виде.

Также ограничим диапазон отображения графика зависимости давления от времени - чтобы не показывать неверные значения получаемых в ходе расчета для малых времен.

# преобразование sympy выражения в lambda функцию с использованием numpy
soln_re_const_p_numpy = sp.lambdify([u, r_d, r_e], soln_sympy.rhs, 
                                    modules = ['numpy', 'scipy'])

# после преобразования lambda функцию можно использовать 
# для проведения расчетов как обычныую функцию

# определим функцию для расчета обратного преобразования Лапласа
# отдельная функция удобно чтобы передавать r_d второй параметр
def soln_re_const_p_inv_lapl_numpy(t_d=1000, r_d=1, r_e=100):
    sol_lap = anaflow.get_lap_inv(soln_re_const_p_numpy,  r_d=r_d, r_e=r_e)
    return sol_lap(t_d).real

# зададим диапазоны изменения параметров для построения графиков
td_arr = np.logspace(1,5,100)
rd_arr = np.logspace(1,2,100)
re=100
# построим графики
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[15,5])
for ri in np.linspace(1,100,10):
    ax1.plot(td_arr, soln_re_const_p_inv_lapl_numpy(td_arr, r_d=ri, r_e=re), label = f'$r_d={ri:.0f}$')
for ti in np.linspace(1e2, 1e4,10):
    ax2.plot(rd_arr, np.squeeze( [soln_re_const_p_inv_lapl_numpy(ti, r_d=ri, r_e=re) for ri in rd_arr ] ), label = f'$t_d={ti:.0f}$')
# наведем немного красоты
ax1.set_xlabel('$t_d$')
ax1.set_ylabel('$p_d$')
ax1.set_ylim(0,5)
ax1.legend()
ax2.set_xlabel('$r_d$')
ax2.set_ylabel('$p_d$')
ax2.legend()
plt.show()

Графики построились, хотя и пришлось подавлять предупреждения об игнорировании мнимой части и для одного из значений расстояния на расчет развалился. Первый вариант построения решения с использованием \(K_0\) выглядит надежнее.

5.7 Частное решение для конечного круглого пласта с замкнутой границей и конечного радиуса скважины.

Частное решение строится за счет поиска параметров \(A\) и \(B\) для уравнения (5.4) так, чтобы они удовлетворяли граничным условиям.

1. Граничное условия на внешнем контуре с замкнутой границей

\[ \left( \frac{d p_D}{d r_d}\right)_{r_d=r_e} = 0 \]

в пространстве Лапласа преобразуется в следующее

\[ \left( \frac{ \tilde{p}_D(r_d, u)}{d r_d}\right)_{r_d=r_e} = 0 \]

2. Граничное условие на скважине

\[ \lim_{r_D \to r_{wD}} \left[ r_D \dfrac{ \partial p_D(r_D, t_D)}{\partial r_D} \right] = -1 \]

в пространстве Лапласа с учетом выражения преобразуется к

\[ \lim_{r_D \to r_{wD}} \left[ r_D \dfrac{ d \tilde{p}_D}{d r_D} \right] = -\dfrac{1}{u} \]

r_e = sp.symbols('r_e', positive=True)
eq1 = sp.Eq(soln2.rhs.diff(r_d).subs(r_d, r_e),0)
eq1

\(\displaystyle - A \sqrt{u} K_{1}\left(r_{e} \sqrt{u}\right) + B \sqrt{u} I_{1}\left(r_{e} \sqrt{u}\right) = 0\)

eq2 = sp.Eq(r_d * soln2.rhs.diff(r_d), -1/u).subs(r_d, 1)
eq2

\(\displaystyle - A \sqrt{u} K_{1}\left(\sqrt{u}\right) + B \sqrt{u} I_{1}\left(\sqrt{u}\right) = - \frac{1}{u}\)

sol_AB = sp.solve([eq1,eq2], [A,B])
display(sol_AB[A])
display(sol_AB[B])

\(\displaystyle \frac{I_{1}\left(r_{e} \sqrt{u}\right)}{- u^{\frac{3}{2}} I_{1}\left(\sqrt{u}\right) K_{1}\left(r_{e} \sqrt{u}\right) + u^{\frac{3}{2}} I_{1}\left(r_{e} \sqrt{u}\right) K_{1}\left(\sqrt{u}\right)}\)

\(\displaystyle \frac{K_{1}\left(r_{e} \sqrt{u}\right)}{- u^{\frac{3}{2}} I_{1}\left(\sqrt{u}\right) K_{1}\left(r_{e} \sqrt{u}\right) + u^{\frac{3}{2}} I_{1}\left(r_{e} \sqrt{u}\right) K_{1}\left(\sqrt{u}\right)}\)

sol_re_close = soln2.subs(A, sol_AB[A]).subs(B, sol_AB[B]).simplify()  
sol_re_close

\(\displaystyle p_{d}{\left(r_{d} \right)} = \frac{I_{0}\left(r_{d} \sqrt{u}\right) K_{1}\left(r_{e} \sqrt{u}\right) + I_{1}\left(r_{e} \sqrt{u}\right) K_{0}\left(r_{d} \sqrt{u}\right)}{u^{\frac{3}{2}} \left(- I_{1}\left(\sqrt{u}\right) K_{1}\left(r_{e} \sqrt{u}\right) + I_{1}\left(r_{e} \sqrt{u}\right) K_{1}\left(\sqrt{u}\right)\right)}\)

Таким образом получено решение для замкнутой круговой границы вокруг скважины \[ \displaystyle \operatorname{p_{d}}{\left(r_{d} \right)} = \frac{I_{0}\left(r_{d} \sqrt{u}\right) K_{1}\left(r_{e} \sqrt{u}\right) + I_{1}\left(r_{e} \sqrt{u}\right) K_{0}\left(r_{d} \sqrt{u}\right)}{u^{\frac{3}{2}} \left(- I_{1}\left(\sqrt{u}\right) K_{1}\left(r_{e} \sqrt{u}\right) + I_{1}\left(r_{e} \sqrt{u}\right) K_{1}\left(\sqrt{u}\right)\right)} \]

Ипользуя подход показанный выше можно построить графики зависимости безразмерного давления от времени и расстояния для этого решения с использованием numpy и scipy функций

# преобразование sympy выражения в lambda функцию с использованием numpy
soln_re_close_numpy = sp.lambdify([u, r_d, r_e], sol_re_close.rhs, 
                                  modules = ['numpy', 'scipy'])

# после преобразования lambda функцию можно использовать 
# для проведения расчетов как обычныую функцию

# определим функцию для расчета обратного преобразования Лапласа
# отдельная функция удобно чтобы передавать r_d второй параметр
def soln_re_close_inv_lapl_numpy(t_d=1000, r_d=1, r_e=100):
    sol_lap = anaflow.get_lap_inv(soln_re_close_numpy,  r_d=r_d, r_e=r_e)
    return sol_lap(t_d)

# зададим диапазоны изменения параметров для построения графиков
td_arr = np.logspace(1,5,100)
rd_arr = np.logspace(1,2,100)
re=100
# построим графики
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=[15,5])
for ri in np.linspace(1,100,10):
    ax1.plot(td_arr, soln_re_close_inv_lapl_numpy(td_arr, r_d=ri, r_e=re), label = f'$r_d={ri:.0f}$')
for ti in np.linspace(1e2, 1e4,10):
    ax2.plot(rd_arr, np.squeeze( [soln_re_close_inv_lapl_numpy(ti, r_d=ri, r_e=re) for ri in rd_arr ] ), label = f'$t_d={ti:.0f}$')
# наведем немного красоты
ax1.set_xlabel('$t_d$')
ax1.set_ylabel('$p_d$')
ax1.legend()
ax2.set_xlabel('$r_d$')
ax2.set_ylabel('$p_d$')
ax2.legend()
plt.show()